Modelling with R

This workflow is based on using the R package kwb.hydrus1d, which is compatible with the last official release of HYDRUS1D v4.17.0140. The development of this new R package was necessary, due to the fact that already available model wrappers were developed for outdated HYDRUS1D versions (e.g. python package phydrus requires HYDRUS1D v4.08) and/or are not well maintained (e.g. open issues/pull requests for R package hydrusR).

The workflow below describes all the steps required for:

Input data preparation
Modifying (i.e. atmospheric input data defined in file ATMOSPHERE.in) an already pre-prepared HYDRUS1D model (using the HYDRUS1D GUI)
Running it from within R by using the command line (cmd) and
Importing and analysing the HYDRUS1D results within R.

Install R Package

# Enable this universe
options(repos = c(
  kwbr = 'https://kwb-r.r-universe.dev',
  CRAN = 'https://cloud.r-project.org'))
# Install R package
install.packages('flextreat.hydrus1d')

Define Paths

paths_list <- list(
  extdata = system.file("extdata", package = "flextreat.hydrus1d"),
  exe_dir = "<extdata>/model",
  model_name = "test",
  model_dir = "<exe_dir>/<model_name>",
  atmosphere = "<model_dir>/ATMOSPH.IN",
  a_level = "<model_dir>/A_LEVEL.out",
  t_level = "<model_dir>/T_LEVEL.out",
  runinf = "<model_dir>/Run_Inf.out",
  solute_id = 1, 
  solute = "<model_dir>/solute<solute_id>.out",
  soil_data = "<extdata>/input-data/soil/soil_geolog.csv"
)

paths <- kwb.utils::resolve(paths_list)

Input Data

Soil

Soil data is derived from depth-dependent grain-size analysis of soil samples taken in Braunschweig. The following required input parameters for the van Genuchten model used in HYDRUS1D were derived by using the pedotransfer function ROSETTA-API version 3, which was developed by Zhang and Schaap, 2017.


library(flextreat.hydrus1d)

soil_dat <- readr::read_csv(paths$soil_data, show_col_types = FALSE) %>% 
  dplyr::mutate(layer_id = dplyr::if_else(.data$tiefe_cm_ende < 60, 
                                       1, 
                                       2),
                thickness_cm = .data$tiefe_cm_ende - .data$tiefe_cm_start,
                sand = .data$S_prozent + .data$G_prozent, 
                silt = .data$U_prozent, 
                clay = round(100 - .data$sand - .data$silt,
                                     1)) %>% 
  dplyr::mutate(clay = dplyr::if_else(.data$clay < 0,
                                      0, 
                                      .data$clay)) %>% 
  soilDB:::ROSETTA(vars = c("sand", "silt", "clay"),
                   v = "3") %>% 
  dplyr::mutate(alpha = 10^.data$alpha,
                npar = 10^.data$npar, 
                ksat = 10^.data$ksat)


knitr::kable(soil_dat)

tiefe_cm	tiefe_cm_start	tiefe_cm_ende	gluehverlust_prozent	U_prozent	S_prozent	G_prozent	kf_beyer	layer_id	thickness_cm	sand	silt	clay	theta_r	theta_s	alpha	npar	ksat	.rosetta.model	.rosetta.version
0 - 30	0	30	2.7	6.4	93.1	0.4	9.5e-05	1	30	93.5	6.4	0.1	0.0481310	0.3671023	0.0334257	3.014652	716.3720	2	3
30 - 55	30	55	1.3	6.3	93.3	0.4	1.3e-04	1	25	93.7	6.3	0.0	0.0480175	0.3669818	0.0335768	3.053357	739.3399	2	3
55 - 80	55	80	0.6	1.7	97.3	1.1	1.9e-04	2	25	98.4	1.7	0.0	0.0485529	0.3626158	0.0355315	3.951321	1296.6154	2	3
80 - 150	80	150	0.2	1.3	98.6	0.1	2.0e-04	2	70	98.7	1.3	0.0	0.0486320	0.3621635	0.0356286	4.021821	1343.1455	2	3
150 - 190	150	190	0.2	1.0	98.7	0.3	2.2e-04	2	40	99.0	1.0	0.0	0.0486788	0.3618558	0.0357316	4.088764	1388.8655	2	3
190 - 210	190	210	0.3	2.5	97.2	0.3	1.7e-04	2	20	97.5	2.5	0.0	0.0484589	0.3633709	0.0351930	3.764078	1171.5976	2	3

Due to similar soil characteristics, only two different layers (column layer_id) were defined: - layer 1: ranging from 0 cm - 55 cm depth - layer 2: ranging from 55 cm - 210 cm depth

The spatially aggregation for each layer (geometric mean) of the input parameters for the van Genuchten model is performed with the code below and shown in the subsequent table.

soil_dat_per_layer <- soil_dat %>% 
  dplyr::group_by(.data$layer_id) %>% 
  dplyr::summarise(layer_thickness = max(.data$tiefe_cm_ende) - min(.data$tiefe_cm_start),
                   theta_r = sum(.data$theta_r*.data$thickness_cm)/.data$layer_thickness,
                   theta_s = sum(.data$theta_s*.data$thickness_cm)/.data$layer_thickness,
                   alpha = sum(.data$alpha*.data$thickness_cm)/.data$layer_thickness,
                   npar = sum(.data$npar*.data$thickness_cm)/.data$layer_thickness,
                   ksat = sum(.data$ksat*.data$thickness_cm)/.data$layer_thickness)


knitr::kable(soil_dat_per_layer)

layer_id	layer_thickness	theta_r	theta_s	alpha	npar	ksat
1	55	0.0480794	0.3670475	0.0334944	3.032245	726.812
2	155	0.0486090	0.3623128	0.0355833	3.994469	1325.304

Atmospheric Boundary Conditions

In total three different atmospheric input boundary conditions (precipitation, potential evaporation and irrigation) are used within the HYDRUS1D model and described in the following subchapters in more detail.

Rainfall

Rainfall is based on historical hourly raw data downloaded from DWD open-data ftp server for station_id = 662 (i.e. Braunschweig). Data is aggregated within R to daily sums by using the code below:

install.packages("rdwd")
library(dplyr)
rdwd::updateRdwd()

rdwd::findID("Braunschweig")
rdwd::selectDWD(name = "Braunschweig", res = "daily")

url_bs_rain <- rdwd::selectDWD(name = "Braunschweig",
                               res = "hourly",
                               var = "precipitation",
                               per = "historical" )

bs_rain <- rdwd::dataDWD(url_bs_rain)

precipitation_hourly <- rdwd::dataDWD(url_bs_rain) %>%
  dplyr::select(.data$MESS_DATUM, .data$R1) %>%
  dplyr::rename("datetime" = "MESS_DATUM",
                "precipitation_mm" = "R1")

usethis::use_data(precipitation_hourly)


precipitation_daily <- precipitation_hourly %>%
  dplyr::mutate("date" = as.Date(datetime)) %>%
  dplyr::group_by(date) %>%
  dplyr::summarise(rain_mm = sum(precipitation_mm))

usethis::use_data(precipitation_daily)

DT::datatable(head(flextreat.hydrus1d::precipitation_daily))

Potential Evaporation

Is based on daily potential evaporation grids (1km x 1km) downloaded from DWD open-data server, where all grids (i.e. 46) within the Abwasserverregnungsgebiet.shp were selected and spatially aggregated (mean, sd, min, max) by using the code below:

shape_file <- system.file("extdata/input-data/gis/Abwasserverregnungsgebiet.shp",
                          package = "flextreat.hydrus1d")

# Only data of full months can currently be read!
evapo_p <- kwb.dwd::read_daily_data_over_shape(
  file = shape_file,
  variable = "evapo_p",
  from = "201701",
  to = "202012"
)

usethis::use_data(evapo_p)

DT::datatable(head(flextreat.hydrus1d::evapo_p))

Irrigation

Monthly irrigation volumes were provided by Abwasserverband Braunschweig for the time period 2017-01-01 - 2020-12-31 as csv file and separates between two sources (groundwater and clearwater). Data preparation was carried out with the code below:

irrigation_file <- system.file("extdata/input-data/Beregnungsmengen_AVB.csv",
                          package = "flextreat.hydrus1d")


irrigation_area <- rgdal::readOGR(dsn = shape_file)


irrigation <- read.csv2(irrigation_file) %>%
  dplyr::select(- .data$Monat) %>%
  dplyr::rename(irrigation_m3 = .data$Menge_m3,
                source = .data$Typ,
                month = .data$Monat_num,
                year = .data$Jahr) %>%
  dplyr::mutate(date_start = as.Date(sprintf("%d-%02d-01",
                               .data$year,
                               .data$month)),
                days_in_month = as.numeric(lubridate::days_in_month(.data$date_start)),
                date_end =  as.Date(sprintf("%d-%02d-%02d",
                                            .data$year,
                                            .data$month,
                                            .data$days_in_month)),
                source = kwb.utils::multiSubstitute(.data$source,
                                                    replacements = list("Grundwasser" = "groundwater.mmPerDay",
                                                                        "Klarwasser" = "clearwater.mmPerDay")),
                irrigation_cbmPerDay = .data$irrigation_m3/.data$days_in_month,
                irrigation_area_sqm = irrigation_area$area,
                irrigation_mmPerDay = 1000*irrigation_cbmPerDay/irrigation_area_sqm) %>%
  dplyr::select(.data$year,
                .data$month,
                .data$days_in_month,
                .data$date_start,
                .data$date_end,
                .data$source,
                .data$irrigation_mmPerDay,
                .data$irrigation_area_sqm) %>%
  tidyr::pivot_wider(names_from = .data$source,
                     values_from = .data$irrigation_mmPerDay)


usethis::use_data(irrigation)

DT::datatable(head(flextreat.hydrus1d::irrigation))

HYDRUS-1D

Modify Input Files

All model input files were initially setup by using the HYDRUS1D-GUI but the following below were modified manually with the

SELECTOR.in

Soil input data were entered manually in SELECTOR.in for the two layers defined the second table in Chapter: Input Data - Soil.

ATMOSPHERE.in

Based on the atmospheric input data (see Chapter: Atmospheric Boundary Conditions) the ATMOSPHERE.in file for HYDRUS1D is prepared with the code below and starts with the hydrological summer half year (assumption: soil is fully wetted at end of winter half year):


atm <- flextreat.hydrus1d::prepare_atmosphere_data()
atm_selected <- flextreat.hydrus1d::select_hydrologic_years(atm)


atm_prep <- flextreat.hydrus1d::prepare_atmosphere(
  atm_selected, 
  conc_irrig_clearwater = 1, # use as "clearwater" tracer 
  conc_irrig_groundwater = 0, 
  conc_rain = 0)

atmos <- kwb.hydrus1d::write_atmosphere(
  atm = atm_prep,
  MaxAL = nrow(atm_prep),
  round_digits = 6 # increase precision for "clearwater" tracer!
  )

writeLines(atmos, paths$atmosphere)

The selected time period covers 1278 days (2017-05-01 - 2020-10-30), i.e. it covers 7 hydrological half years.

DT::datatable(atm_selected)

These time-series were converted with the function flextreat.hydrus1d::prepare_atmosphere() into the data format required by HYDRUS1D. Due to the fact, that irrigation rates (i.e. sum of clearwater.mmPerDay and groundwater.mmPerDay) cannot be entered separately as input column within HYDRUS1D, these were simply added to th prec (i.e. precipitation) column. The whole time series defined in ATMOSPHERE.in is shown below:

DT::datatable(atm_prep)

Run Model

Finally the model is run automatically by using the following code:

exe_path <- kwb.hydrus1d::check_hydrus_exe(dir = paths$exe_dir,
                                           skip_preinstalled = TRUE)
#> Checking if download of HYDRUS1D executable v4.17.0140 from 'https://github.com/mrustl/hydrus1d/archive/refs/tags/v4.17.0140.zip' was successful ... ok. (0.00s)
kwb.hydrus1d:::run_model(exe_path = exe_path,
                         model_path = paths$model_dir)
#> Warning in shell(cmd = sprintf("cd %s && %s", fs::path_abs(target_dir), : 'cd
#> D:/a/_temp/Library/flextreat.hydrus1d/extdata/model && H1D_CALC.exe' execution
#> failed with error code 24

Read Results


runinf <- kwb.hydrus1d::read_runinf(paths$runinf)

summary(runinf)
#>     t_level           time               dt                itr_w       
#>  Min.   :    1   Min.   :   0.05   Min.   :0.0001065   Min.   : 2.000  
#>  1st Qu.: 5954   1st Qu.: 285.35   1st Qu.:0.0220900   1st Qu.: 2.000  
#>  Median :11908   Median : 649.70   Median :0.0495000   Median : 2.000  
#>  Mean   :11908   Mean   : 624.94   Mean   :0.0537056   Mean   : 2.483  
#>  3rd Qu.:17862   3rd Qu.: 957.35   3rd Qu.:0.0944970   3rd Qu.: 2.000  
#>  Max.   :23815   Max.   :1279.00   Max.   :0.1000000   Max.   :20.000  
#>      itr_c       it_cum          kod_t               kod_b    converg       
#>  Min.   :1   Min.   :    4   Min.   :-4.000000   Min.   :-5   Mode:logical  
#>  1st Qu.:1   1st Qu.:19307   1st Qu.:-4.000000   1st Qu.:-5   TRUE:23815    
#>  Median :1   Median :36910   Median : 4.000000   Median :-5                 
#>  Mean   :1   Mean   :37823   Mean   : 0.004199   Mean   :-5                 
#>  3rd Qu.:1   3rd Qu.:56331   3rd Qu.: 4.000000   3rd Qu.:-5                 
#>  Max.   :1   Max.   :75157   Max.   : 4.000000   Max.   :-5                 
#>      peclet       courant      
#>  Min.   :0.1   Min.   :0.0000  
#>  1st Qu.:0.1   1st Qu.:0.0150  
#>  Median :0.1   Median :0.0400  
#>  Mean   :0.1   Mean   :0.0829  
#>  3rd Qu.:0.1   3rd Qu.:0.1170  
#>  Max.   :0.1   Max.   :0.9960

a_level <- kwb.hydrus1d::read_alevel(paths$a_level)
a_level
#> # A tibble: 1,279 × 10
#>     time sum_r_top sum_r…¹ sum_v…² sum_v…³ sum_v_…⁴   h_top h_root h_bot a_level
#>  * <dbl>     <dbl>   <dbl>   <dbl>   <dbl>    <dbl>   <dbl>  <dbl> <dbl>   <dbl>
#>  1     1    0.246        0  0.0672       0 -0.00173 -1   e5      0 -116.       1
#>  2     2   -0.213        0 -0.392        0 -0.00346 -8.30e1      0 -116.       2
#>  3     3   -0.0375       0 -0.216        0 -0.00520 -1.67e3      0 -116.       3
#>  4     4   -0.162        0 -0.341        0 -0.00693 -1.03e2      0 -116.       4
#>  5     5   -0.157        0 -0.335        0 -0.00866 -1.23e2      0 -116.       5
#>  6     6   -0.0311       0 -0.221        0 -0.0104  -1   e5      0 -116.       6
#>  7     7    0.188        0 -0.187        0 -0.0121  -1   e5      0 -116.       7
#>  8     8    0.232        0 -0.167        0 -0.0139  -1   e5      0 -116.       8
#>  9     9    0.408        0 -0.151        0 -0.0156  -1   e5      0 -116.       9
#> 10    10    0.583        0 -0.139        0 -0.0173  -1   e5      0 -116.      10
#> # … with 1,269 more rows, and abbreviated variable names ¹sum_r_root,
#> #   ²sum_v_top, ³sum_v_root, ⁴sum_v_bot

t_level <- kwb.hydrus1d::read_tlevel(paths$t_level)
t_level
#> # A tibble: 23,815 × 22
#>     time r_top r_root  v_top v_root    v_bot sum_r_top sum_r_r…¹ sum_v…² sum_v…³
#>  * <dbl> <dbl>  <dbl>  <dbl>  <dbl>    <dbl>     <dbl>     <dbl>   <dbl>   <dbl>
#>  1 0.05  0.246      0 0.246       0 -0.00173    0.0123         0  0.0123       0
#>  2 0.1   0.246      0 0.246       0 -0.00173    0.0246         0  0.0246       0
#>  3 0.15  0.246      0 0.246       0 -0.00173    0.0369         0  0.0369       0
#>  4 0.168 0.246      0 0.246       0 -0.00173    0.0414         0  0.0414       0
#>  5 0.188 0.246      0 0.0664      0 -0.00173    0.0464         0  0.0427       0
#>  6 0.203 0.246      0 0.0433      0 -0.00173    0.0498         0  0.0433       0
#>  7 0.218 0.246      0 0.0562      0 -0.00173    0.0536         0  0.0442       0
#>  8 0.235 0.246      0 0.0596      0 -0.00173    0.0578         0  0.0452       0
#>  9 0.254 0.246      0 0.0575      0 -0.00173    0.0625         0  0.0463       0
#> 10 0.275 0.246      0 0.0494      0 -0.00173    0.0675         0  0.0473       0
#> # … with 23,805 more rows, 12 more variables: sum_v_bot <dbl>, h_top <dbl>,
#> #   h_root <dbl>, h_bot <dbl>, run_off <dbl>, sum_run_off <dbl>, volume <dbl>,
#> #   sum_infil <dbl>, sum_evap <dbl>, t_level <dbl>, cum_w_trans <dbl>,
#> #   snow_layer <dbl>, and abbreviated variable names ¹sum_r_root, ²sum_v_top,
#> #   ³sum_v_root

solute <- kwb.hydrus1d::read_solute(paths$solute)
solute
#> # A tibble: 23,815 × 19
#>     time cv_top cv_bot sum_cv…¹ sum_c…² cv_ch0 cv_ch1 c_top c_root c_bot cv_root
#>  * <dbl>  <dbl>  <dbl>    <dbl>   <dbl>  <dbl>  <dbl> <dbl>  <dbl> <dbl>   <dbl>
#>  1 0.05   1.57       0   0.0783       0      0      0  1.00      0     0       0
#>  2 0.1    0.779      0   0.117        0      0      0  1.00      0     0       0
#>  3 0.15   0.558      0   0.145        0      0      0  1.00      0     0       0
#>  4 0.168  0.444      0   0.153        0      0      0  1.00      0     0       0
#>  5 0.188  0.298      0   0.159        0      0      0  1.00      0     0       0
#>  6 0.203  0.143      0   0.161        0      0      0  1.00      0     0       0
#>  7 0.218  0.129      0   0.163        0      0      0  1.00      0     0       0
#>  8 0.235  0.138      0   0.166        0      0      0  1.00      0     0       0
#>  9 0.254  0.134      0   0.168        0      0      0  1.00      0     0       0
#> 10 0.275  0.120      0   0.171        0      0      0  1.00      0     0       0
#> # … with 23,805 more rows, 8 more variables: sum_cv_root <dbl>,
#> #   sum_cv_n_eql <dbl>, t_level <dbl>, c_gwl <dbl>, c_run_off <dbl>,
#> #   sum_c_run_off <dbl>, cv_i <dbl>, sum_cv_i <dbl>, and abbreviated variable
#> #   names ¹sum_cv_top, ²sum_cv_bot

Plot


gg <- solute %>%
  dplyr::select(.data$time, .data$c_top, .data$c_bot) %>%
  dplyr::mutate(date = as.POSIXct("2017-05-01", tz = "UTC") + .data$time*24*3600) %>%
  tidyr::pivot_longer(cols = c("c_top", "c_bot"),
                      names_to = "parameter",
                      values_to = "value") %>%
  ggplot2::ggplot(mapping = ggplot2::aes(x = .data$date,
                                         y = .data$value,
                                         col = .data$parameter )) +
  ggplot2::geom_point() +
  ggplot2::geom_line() +
  ggplot2::ylim(c(0,1)) +
  ggplot2::theme_bw()
gg

plotly::ggplotly(gg)

Michael Rustler